*************************************************************************************************
*																								*
*						Flexible Wages, Bargaining, and The Gender Gap							*
*						Barbara Biasi and Heather Sarsons										*
*						Abraham and Sun (2020) test												*
*																								*
*************************************************************************************************

/*
clear all
set more off
set matsize 11000
set maxvar 32000


global user = 1 // 1 = Heather, 2 = Barbara

if $user == 1 {
cd "/Users/sarsons/Dropbox/wisconsin_women/data"
global tab = "/Users/sarsons/Dropbox/wisconsin_women/draft/tables"
global out = "/Users/sarsons/Dropbox/wisconsin_women/out"
global RA = "/Users/sarsons/Dropbox/wisconsin_women/do/calvin"
global data "/Users/sarsons/Dropbox/wisconsin_women/data"
}
if $user == 2 {
cd "~/Dropbox/Research/wisconsin_women/data"
global tab = "~/Dropbox/Research/wisconsin_women/draft/tables"
global out = "~/Dropbox/Research/wisconsin_women/out"
global RA = "~/Dropbox/Research/wisconsin_women/do/calvin"
}

* Yale colors
global yaleblue = "0 53 107"
global ylb = "40 109 192"
global yo = "189 83 25"
global ylight = "217 233 242"
global graph = "xlabel(,grid glp(dot) glc(gs10) tlc(gs10)) ylabel(,grid glp(dot) glc(gs10) tlc(gs10)) plotregion(lp(blank))"
global graphbar = "ylabel(,grid glp(dot) glc(gs10) tlc(gs10)) plotregion(lp(blank))"


* Load data
use teachers_panel_va.dta, clear

* Sample
keep if year >= 2007 & year <= 2016 

* new agreement data
drop expire extension
sort district_code
merge m:1 district_code using new_agreements.dta
drop if _m == 2
drop _m

* the following change sets exp = 2011 even to districts w/first handbook in 2012.

replace expire = 2011 if doubt == "possibly 2011"
replace extension = 0 if doubt == "possibly 2011"
replace expire = 2011 if doubt == "possibly exp = 2011, no ext"
replace extension = 0 if doubt == "possibly exp = 2011, no ext"
replace expire = 2013 if doubt == "possibly exp = 2013, no ext"
replace extension = 0 if doubt == "possibly exp = 2013, no ext"


*
drop doubt
gen Extension = extension
replace extension = expire if extension == 0



* Set extension to expiry date if no extension
* Set extension to expiry date if no extension
replace extension = expire if extension == .

* Merge admin info
preserve
use administrators.dta, clear
keep if year > 2006 & year <= 2016
keep if position == 51 // & year == 2011
collapse male_princ = female, by(schoolcode district_code year)
replace male_princ = round(male)
replace male_princ = 1 - male_princ
sort schoolcode district_code year
save temp.dta, replace
restore
capt rename school_code schoolcode
sort schoolcode district_code year
merge m:1 schoolcode district_code year using temp.dta
drop _m
rm temp.dta
preserve
use administrators.dta, clear
keep if year > 2006 & year <= 2016
keep if position == 5 // & year == 2011
collapse male_super = female, by(district_code year)
replace male_super = round(male)
replace male_super = 1 - male_super
sort district_code year
save temp.dta, replace
restore
sort district_code year
merge m:1 district_code year using temp.dta
drop _m
rm temp.dta

* principal during exp year
gen P = male_princ if year == extension
bysort id: egen maleprinc = max(P)
drop P
* super during exp year
gen P = male_super if year == extension
bysort id: egen malesuper = max(P)
drop P

* mean gender principal in the years before extension during exp year
gen P = male_princ if year <= extension
bysort id: egen maleprinc_pre = mean(P)
replace maleprinc_pre = round(maleprinc_pre)
drop P
* mean gender super in the years before extension during exp year
gen P = male_super if year <= extension
bysort id: egen malesuper_pre = mean(P)
replace malesuper_pre = round(malesuper_pre)
drop P

* Interaction

g male_super_princ = (malesuper_pre==1&maleprinc_pre==1) 
g fem_super_princ = (malesuper_pre==0&maleprinc_pre==0)
g fems_malep = (malesuper_pre==0&maleprinc_pre==1)
g males_femp = (malesuper_pre==1&maleprinc_pre==0)

* moves by CZ			
	merge m:1 district_code using district_county_cz.dta
		drop if _m==2
		drop _m
		
	sort id year
	by id: g moved_cz = 1 if cz[_n]!=cz[_n-1]
	replace moved_cz=0 if moved_cz==.	
	by id: replace moved_cz=. if _n==1	
	
		* Moves within CZ
			sort id year
			by id: g moves_within_cz = 1 if (cz==cz[_n-1]) & mover_d == 1
			replace moves_within_cz=0 if moves_within_cz==.	
			by id: replace moves_within_cz=. if _n==1
			
replace moves_within_cz=0 if mover_d==0	
gen moves_x_cz = mover_d == 1 & moves_within_cz == 0
replace moves_x_cz = . if moves_within_cz == . | mover_d == .

	
* ever moves
	bys id: egen ever_moves = max(mover_d)	
	
tempvar x
gen `x' = year if mover_d == 1
bysort id: egen mover_year = max(`x')
gen timetomove = year - mover_year
local x = 1
qui sum timetomove
forvalues n = `r(min)'/ `r(max)' {
gen move`x' = timetomove == `n'
label var move`x' "`n'"
local x = `x' + 1
}	


*--------->  Variables we need

gen logsalary = log(salary_n)

gen postexp = year > expire
replace postexp = . if expire == .

gen postext = year > ext
replace postext = . if ext == .

gen masterup = master == 1 | phd == 1


* controls: grade
gen kinder 	= highgrade == "K3" | highgrade == "K4" | highgrade == "KG" | highgrade == "PK"
gen elem 	= highgrade == "01" | highgrade == "02" | highgrade == "03" | highgrade == "04" | highgrade == "05" | highgrade == "06"
gen middle 	= highgrade == "07" | highgrade == "08"
gen high 	= highgrade == "09" | highgrade == "10" | highgrade == "11" | highgrade == "12"


* dummy for teachers in top and bottom quartile of value-added distribution
qui sum va_s, det
gen highva = 1 if va_s > r(p50)
replace highva = 0 if va_s <= r(p50)
replace highva = . if va_s == .

* absorbed controls
global exp = "i.district_code i.district_code#i.postexp i.totalexp i.totalexp#i.postexp i.master i.master#i.postexp i.phd i.phd#i.postexp i.high i.high#i.postexp i.math i.math#i.postexp i.year##i.expire"
global exp_pp = "i.district_code i.pp i.district_code#i.postexp i.totalexp i.pp#i.totalexp i.totalexp#i.postexp i.pp#i.totalexp#i.postexp i.masterup i.pp#i.masterup i.masterup#i.postexp i.pp#i.masterup#i.postexp i.high i.pp#i.high i.high#i.postexp i.pp#i.high#i.postexp i.year i.pp#i.year i.year#i.expire i.pp#i.year#i.expire"
global ext = "i.district_code i.district_code#i.postext i.totalexp i.totalexp#i.postext i.master i.master#i.postext i.phd i.phd#i.postext i.high i.high#i.postext i.math i.math#i.postext"
*--> time variables

* expiration year
gen time = year - expire
qui tab time 
forvalues n = 1/ `r(r)' {
local z = `n' - 7
gen fem_`n' = female * (time == `z')
replace fem_`n' = . if time == .
label var fem_`n' "`z'"
}

* generate trick variables for table labeling

gen Ed = 0
gen Edexp = 0
gen Edext = 0
gen Yr = 0

global tag_exp = "Ed Edexp Yr"
global tag_ext = "Ed Edext Yr"

label var Ed "District, Educ, Exper"
label var Edexp "District, Educ, Exper * post-expir"
label var Edext "District, Educ, Exper * post-extens"
label var Yr "Year"
*/

*-------------------> ABRAHAM-SUN SET UP <-------------------*

g reltime = year-extension
g l_5 = (reltime==-5)
g l_4 = (reltime==-4)
g l_3 = (reltime==-3)
g l_2 = (reltime==-2)
g l_1 = (reltime==-1)
g l0 = (reltime==0)
g l1 = (reltime==1)
g l2 = (reltime==2)
g l3 = (reltime==3)
g l4 = (reltime==4)
g l5 = (reltime==5)

*bys id: egen treated = min(extension)
g t = extension if year==2011
bys id: egen treated = max(t)
drop t
g t1 = treated==2011
g t2 = treated==2012
g t3 = treated==2013
g t4 = treated==2014
g t5 = treated==2016

/* UNCOMMENT TO CALCULATE EVENT STUDY WEIGHTS
eventstudyweights l_5 l_4 l_3 l_2 l_1 l0 l1 l2 l3 l4 l5, controls(i.district_code i.year) cohort(treated) rel_time(reltime) saveweights("weights")

preserve
	insheet using weights.csv, clear

	keep l_2 treated reltime
	reshape wide l_2, i(reltime) j(treated)

	drop if reltime<-5
	graph twoway line l_2* reltime, legend(order(1 "2011" 2 "2012" 3 "2013" 4 "2014") rows(1)) ///
		xla(-5(1)5) graphreg(fc(white)) yla(-1(0.5)1) xtitle("Relative Time") ytitle("Weights")
	graph export $out/weights_AS2020.png, replace
	egen w_sum = rowtotal(l_2*)
restore	
*/

* Generate counts for each cohort
count if treated==2011 & year==2011
global c_2011 = r(N)
count if treated==2012 & year==2011
global c_2012 = r(N)
count if treated==2013 & year==2011
global c_2013 = r(N)
count if treated==2014 & year==2011
global c_2014 = r(N)
count if treated==2016 & year==2011
global c_2016 = r(N)

*qui tab reltime
gen timex = year - ext
qui tab timex
forvalues n = 1/ `r(r)' {
local z = `n' - 8
gen femx_`n' = female * (reltime == `z')
replace femx_`n' = . if reltime == .
label var femx_`n' "`z'"
}
egen Timex = group(reltime)
g zero=0
label var zero "0"

save $data/AS_setup.dta, replace
	
*----------------> BY AGE

	use $data/AS_setup.dta, clear
	gen age = year - birth
	g young = 1 if (age<=34)
	replace young = 0 if(age>=52&age<=65)

	keep if young==1

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

	reghdfe logsalary $rhs_rel_year_i, a($exp i.year##i.extension) resid

	mat b = e(b)
	mat V = vecdiag(e(V))
	mat b_2011 = (b[.,1..10])
	mat V_2011 = (V[.,1..10])	
	mat b_2012 = (b[.,11..20])
	mat V_2012 = (V[.,11..20])
	mat b_2013 = (b[.,21..30])
	mat V_2013 = (V[.,21..30])
	mat b_2014 = (b[.,31..40])
	mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2
		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014

	// output results to excel spreadsheet
	clear
	global n_obs = 11
	set obs $n_obs
	gen y = _n-5 // generate relative year

	* two-way FE estimates
	*mat b2way_hrs = b2way_hrs[1..4,.]\(0,0)\b2way_hrs[5..9,.] // fill in zeros for evt_time_4, which are normalized to zero
	svmat b2way_hrs 

	gen ci1_u = sqrt(b2way_hrs2)*1.96+b2way_hrs1
	gen ci1_l = -sqrt(b2way_hrs2)*1.96+b2way_hrs1
	gen sd1 = sqrt(b2way_hrs2)

	* IW estimates
	svmat biw_hrs

	gen ci2_u = sqrt(biw_hrs2)*1.96+biw_hrs1
	gen ci2_l = -sqrt(biw_hrs2)*1.96+biw_hrs1
	gen sd2 = sqrt(biw_hrs2)
	* CATTs estimates are biw_hrs3 - biw_hrs7
	gen sd2011 = sqrt(biw_hrs4)
	gen sd2012 = sqrt(biw_hrs6)
	gen sd2013 = sqrt(biw_hrs8)
	gen sd2014 = sqrt(biw_hrs10)
	order y b2way_hrs1 biw_hrs1 biw_hrs3 biw_hrs5 biw_hrs7 biw_hrs9 sd1 sd2 sd2011 sd2012 sd2013 sd2014
		
	export delimited using $data/weighted_age.csv, replace
		
		
		
	* Older teachers
		
	use $data/AS_setup.dta, clear
	gen age = year - birth
	g young = 1 if (age<=34)
	replace young = 0 if(age>=52&age<=65)

	keep if young==0

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs_old = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2
		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs_old = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs_old 

		gen ci1_u_old = sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen ci1_l_old = -sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen sd1_old = sqrt(b2way_hrs_old2)

		* IW estimates
		svmat biw_hrs_old

		gen ci2_u_old = sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen ci2_l_old = -sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen sd2_old = sqrt(biw_hrs_old2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011_old = sqrt(biw_hrs_old4)
		gen sd2012_old = sqrt(biw_hrs_old6)
		gen sd2013_old = sqrt(biw_hrs_old8)
		gen sd2014_old = sqrt(biw_hrs_old10)
		order y b2way_hrs_old1 biw_hrs_old1 biw_hrs_old3 biw_hrs_old5 biw_hrs_old7 biw_hrs_old9 sd1_old sd2_old sd2011_old sd2012_old sd2013_old sd2014_old
		
		export delimited using $data/weighted_age2.csv, replace
		
	
		appendfile $data/weighted_age.csv $data/weighted_age2.csv

*----------------> BY EXPERIENCE

	use $data/AS_setup.dta, clear
	qui sum totalexp, d
	g low_exp = 1 if (totalexp<=`r(p25)')
	replace low_exp = 0 if(totalexp>=`r(p75)' & totalexp!=.)

	keep if low_exp==1

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2
	
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs 

		gen ci1_u = sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen ci1_l = -sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen sd1 = sqrt(b2way_hrs2)

		* IW estimates
		svmat biw_hrs

		gen ci2_u = sqrt(biw_hrs2)*1.96+biw_hrs1
		gen ci2_l = -sqrt(biw_hrs2)*1.96+biw_hrs1
		gen sd2 = sqrt(biw_hrs2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011 = sqrt(biw_hrs4)
		gen sd2012 = sqrt(biw_hrs6)
		gen sd2013 = sqrt(biw_hrs8)
		gen sd2014 = sqrt(biw_hrs10)
		order y b2way_hrs1 biw_hrs1 biw_hrs3 biw_hrs5 biw_hrs7 biw_hrs9 sd1 sd2 sd2011 sd2012 sd2013 sd2014
		
		export delimited using $data/weighted_exp1.csv, replace
		
		
		
	* More experience
		
	use $data/AS_setup.dta, clear
	qui sum totalexp, d
	g low_exp = 1 if (totalexp<=`r(p25)')
	replace low_exp = 0 if(totalexp>=`r(p75)' & totalexp!=.)

	keep if low_exp==0

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs_old = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2

		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs_old = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs_old 

		gen ci1_u_old = sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen ci1_l_old = -sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen sd1_old = sqrt(b2way_hrs_old2)

		* IW estimates
		svmat biw_hrs_old

		gen ci2_u_old = sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen ci2_l_old = -sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen sd2_old = sqrt(biw_hrs_old2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011_old = sqrt(biw_hrs_old4)
		gen sd2012_old = sqrt(biw_hrs_old6)
		gen sd2013_old = sqrt(biw_hrs_old8)
		gen sd2014_old = sqrt(biw_hrs_old10)
		order y b2way_hrs_old1 biw_hrs_old1 biw_hrs_old3 biw_hrs_old5 biw_hrs_old7 biw_hrs_old9 sd1_old sd2_old sd2011_old sd2012_old sd2013_old sd2014_old
		
		export delimited using $data/weighted_exp2.csv, replace
			
		appendfile $data/weighted_exp1.csv $data/weighted_exp2.csv
		
	


*----------------> BY PRINCIPAL GENDER

	* Male Principal
	
	use $data/AS_setup.dta, clear

	keep if maleprinc_pre == 1

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2
		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs 

		gen ci1_u = sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen ci1_l = -sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen sd1 = sqrt(b2way_hrs2)

		* IW estimates
		svmat biw_hrs

		gen ci2_u = sqrt(biw_hrs2)*1.96+biw_hrs1
		gen ci2_l = -sqrt(biw_hrs2)*1.96+biw_hrs1
		gen sd2 = sqrt(biw_hrs2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011 = sqrt(biw_hrs4)
		gen sd2012 = sqrt(biw_hrs6)
		gen sd2013 = sqrt(biw_hrs8)
		gen sd2014 = sqrt(biw_hrs10)
		order y b2way_hrs1 biw_hrs1 biw_hrs3 biw_hrs5 biw_hrs7 biw_hrs9 sd1 sd2 sd2011 sd2012 sd2013 sd2014
		
		export delimited using $data/weighted_princ1.csv, replace
		
		
		
	* Female Principal
		
	use $data/AS_setup.dta, clear
	
	keep if maleprinc_pre == 0

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs_old = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
			
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2

		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs_old = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs_old 

		gen ci1_u_old = sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen ci1_l_old = -sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen sd1_old = sqrt(b2way_hrs_old2)

		* IW estimates
		svmat biw_hrs_old

		gen ci2_u_old = sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen ci2_l_old = -sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen sd2_old = sqrt(biw_hrs_old2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011_old = sqrt(biw_hrs_old4)
		gen sd2012_old = sqrt(biw_hrs_old6)
		gen sd2013_old = sqrt(biw_hrs_old8)
		gen sd2014_old = sqrt(biw_hrs_old10)
		order y b2way_hrs_old1 biw_hrs_old1 biw_hrs_old3 biw_hrs_old5 biw_hrs_old7 biw_hrs_old9 sd1_old sd2_old sd2011_old sd2012_old sd2013_old sd2014_old
		
		export delimited using $data/weighted_princ2.csv, replace
			
		appendfile $data/weighted_princ1.csv $data/weighted_princ2.csv
		


*----------------> BY SUPERINTENDENT GENDER

	* Male Super
	
	use $data/AS_setup.dta, clear

	keep if male_super_year == 1

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2
			
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs 

		gen ci1_u = sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen ci1_l = -sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen sd1 = sqrt(b2way_hrs2)

		* IW estimates
		svmat biw_hrs

		gen ci2_u = sqrt(biw_hrs2)*1.96+biw_hrs1
		gen ci2_l = -sqrt(biw_hrs2)*1.96+biw_hrs1
		gen sd2 = sqrt(biw_hrs2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011 = sqrt(biw_hrs4)
		gen sd2012 = sqrt(biw_hrs6)
		gen sd2013 = sqrt(biw_hrs8)
		gen sd2014 = sqrt(biw_hrs10)
		order y b2way_hrs1 biw_hrs1 biw_hrs3 biw_hrs5 biw_hrs7 biw_hrs9 sd1 sd2 sd2011 sd2012 sd2013 sd2014
		
		export delimited using $data/weighted_super1.csv, replace
		
		
		
	* Female Super
		
	use $data/AS_setup.dta, clear
	
	keep if male_super_year == 0

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs_old = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}
		
	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"
	reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 



		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs_old = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs_old 

		gen ci1_u_old = sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen ci1_l_old = -sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen sd1_old = sqrt(b2way_hrs_old2)

		* IW estimates
		svmat biw_hrs_old

		gen ci2_u_old = sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen ci2_l_old = -sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen sd2_old = sqrt(biw_hrs_old2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011_old = sqrt(biw_hrs_old4)
		gen sd2012_old = sqrt(biw_hrs_old6)
		gen sd2013_old = sqrt(biw_hrs_old8)
		gen sd2014_old = sqrt(biw_hrs_old10)
		order y b2way_hrs_old1 biw_hrs_old1 biw_hrs_old3 biw_hrs_old5 biw_hrs_old7 biw_hrs_old9 sd1_old sd2_old sd2011_old sd2012_old sd2013_old sd2014_old
		
		export delimited using $data/weighted_super2.csv, replace
			
		appendfile $data/weighted_super1.csv $data/weighted_super2.csv
	

*----------------> BY DISTRICT TYPE (APPENDIX)

	* No schedule
	
	use $data/AS_setup.dta, clear

	keep if pp == 1

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2		
		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs 

		gen ci1_u = sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen ci1_l = -sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen sd1 = sqrt(b2way_hrs2)

		* IW estimates
		svmat biw_hrs

		gen ci2_u = sqrt(biw_hrs2)*1.96+biw_hrs1
		gen ci2_l = -sqrt(biw_hrs2)*1.96+biw_hrs1
		gen sd2 = sqrt(biw_hrs2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011 = sqrt(biw_hrs4)
		gen sd2012 = sqrt(biw_hrs6)
		gen sd2013 = sqrt(biw_hrs8)
		gen sd2014 = sqrt(biw_hrs10)
		order y b2way_hrs1 biw_hrs1 biw_hrs3 biw_hrs5 biw_hrs7 biw_hrs9 sd1 sd2 sd2011 sd2012 sd2013 sd2014
		
		export delimited using $data/weighted_pp1.csv, replace
		
		
		
	* Female Principal
		
	use $data/AS_setup.dta, clear
	
	keep if pp == 0

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs_old = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2		

		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs_old = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs_old 

		gen ci1_u_old = sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen ci1_l_old = -sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen sd1_old = sqrt(b2way_hrs_old2)

		* IW estimates
		svmat biw_hrs_old

		gen ci2_u_old = sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen ci2_l_old = -sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen sd2_old = sqrt(biw_hrs_old2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011_old = sqrt(biw_hrs_old4)
		gen sd2012_old = sqrt(biw_hrs_old6)
		gen sd2013_old = sqrt(biw_hrs_old8)
		gen sd2014_old = sqrt(biw_hrs_old10)
		order y b2way_hrs_old1 biw_hrs_old1 biw_hrs_old3 biw_hrs_old5 biw_hrs_old7 biw_hrs_old9 sd1_old sd2_old sd2011_old sd2012_old sd2013_old sd2014_old
		
		export delimited using $data/weighted_pp2.csv, replace
			
		appendfile $data/weighted_pp1.csv $data/weighted_pp2.csv
		
		

	
*----------------> MOBILITY GRAPHS

	* Never moves
	
	use $data/AS_setup.dta, clear

	keep if mover_year==.

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2		
		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs 

		gen ci1_u = sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen ci1_l = -sqrt(b2way_hrs2)*1.96+b2way_hrs1
		gen sd1 = sqrt(b2way_hrs2)

		* IW estimates
		svmat biw_hrs

		gen ci2_u = sqrt(biw_hrs2)*1.96+biw_hrs1
		gen ci2_l = -sqrt(biw_hrs2)*1.96+biw_hrs1
		gen sd2 = sqrt(biw_hrs2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011 = sqrt(biw_hrs4)
		gen sd2012 = sqrt(biw_hrs6)
		gen sd2013 = sqrt(biw_hrs8)
		gen sd2014 = sqrt(biw_hrs10)
		order y b2way_hrs1 biw_hrs1 biw_hrs3 biw_hrs5 biw_hrs7 biw_hrs9 sd1 sd2 sd2011 sd2012 sd2013 sd2014
		
		export delimited using $data/weighted_mover1.csv, replace
		
		
		
	* Ever movers
		
	use $data/AS_setup.dta, clear
	
	keep if mover_year != .

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_13 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs_old = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2		
		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs_old = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs_old 

		gen ci1_u_old = sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen ci1_l_old = -sqrt(b2way_hrs_old2)*1.96+b2way_hrs_old1
		gen sd1_old = sqrt(b2way_hrs_old2)

		* IW estimates
		svmat biw_hrs_old

		gen ci2_u_old = sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen ci2_l_old = -sqrt(biw_hrs_old2)*1.96+biw_hrs_old1
		gen sd2_old = sqrt(biw_hrs_old2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011_old = sqrt(biw_hrs_old4)
		gen sd2012_old = sqrt(biw_hrs_old6)
		gen sd2013_old = sqrt(biw_hrs_old8)
		gen sd2014_old = sqrt(biw_hrs_old10)
		order y b2way_hrs_old1 biw_hrs_old1 biw_hrs_old3 biw_hrs_old5 biw_hrs_old7 biw_hrs_old9 sd1_old sd2_old sd2011_old sd2012_old sd2013_old sd2014_old
		
		export delimited using $data/weighted_mover2.csv, replace
		
	* Post movers
		
	use $data/AS_setup.dta, clear
	
	keep if mover_year != .& mover_year > extension

	reghdfe logsalary femx_5-femx_9 zero femx_11-femx_15 female if reltime > -7, a($exp i.extension##i.year) vce(cluster district_code)
		
		// store results
		mat b = e(b)
		mat b = b[.,1..11]
		mat V = vecdiag(e(V))
		mat V = V[.,1..11]	
		mat b2way_hrs_m = b', V'

		// exclude the last cohort 
		foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
			replace `vv' = 0 if treated == 2016
		}

		preserve
			foreach vv of varlist femx_5-femx_9 zero femx_11-femx_15 {
				replace `vv' = . if treated == 2016
			}  
			qui reghdfe t1 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res1, resid		
			qui reghdfe t2 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res2, resid		
			qui reghdfe t3 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res3, resid
			qui reghdfe t4 if treated<2016&reltime>-7, a($exp i.year) resid
				predict res4, resid
				
			reg res1 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r1
			reg res2 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r2
			reg res3 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r3
			reg res4 femx_5-femx_9 zero femx_11-femx_15 if treated<2016&reltime>-7
				est store r4
				
			suest r1 r2 r3 r4
			matrix V = e(V)
			matrix list V
			matselrc V Vcov, r(1/10,13/22,25/34,37/46) c(1/10,13/22,25/34,37/46)
		restore

		foreach vv of varlist femx_5-femx_9 femx_11-femx_15 {
			forval i=2011/2016 {
				gen `vv'_`i' = `vv' * (ext==`i')
			}
		}

	global rhs_rel_year_i ""
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2011 femx_6_2011 femx_7_2011 femx_8_2011 femx_9_2011 femx_11_2011 femx_12_2011 femx_13_2011 femx_14_2011 femx_15_2011"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2012 femx_6_2012 femx_7_2012 femx_8_2012 femx_9_2012 femx_11_2012 femx_12_2012 femx_13_2012 femx_14_2012 femx_15_2012"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2013 femx_6_2013 femx_7_2013 femx_8_2013 femx_9_2013 femx_11_2013 femx_12_2013 femx_13_2013 femx_14_2013 femx_15_2013"
	global rhs_rel_year_i "$rhs_rel_year_i femx_5_2014 femx_6_2014 femx_7_2014 femx_8_2014 femx_9_2014 femx_11_2014 femx_12_2014 femx_13_2014 femx_14_2014 femx_15_2014"

		reghdfe logsalary $rhs_rel_year_i if reltime>-7, a($exp i.year##i.extension) resid

		mat b = e(b)
		mat V = vecdiag(e(V))
		mat b_2011 = (b[.,1..10])
		mat V_2011 = (V[.,1..10])	
		mat b_2012 = (b[.,11..20])
		mat V_2012 = (V[.,11..20])
		mat b_2013 = (b[.,21..30])
		mat V_2013 = (V[.,21..30])
		mat b_2014 = (b[.,31..40])
		mat V_2014 = (V[.,31..40])

		mat biw_long_m = b 


		* Coef on event time 5 for cohorts 2012, 2013, 2014
		matselrc biw_long tempb, c(11,21,31)
		matselrc Vcov tempVcov, r(11,21,31) c(11,21,31)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2012/($c_2012+$c_2013+$c_2014)*femx_5_2012 + $c_2013/($c_2012+$c_2013+$c_2014)*femx_5_2013 + $c_2014/($c_2012+$c_2013+$c_2014)*femx_5_2014
		matrix b = r(estimate)
		matrix V = r(se)^2 + temp
		
		* Coef on event time 6 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(2,12,22,32)
		matselrc Vcov tempVcov, r(2,12,22,32) c(2,12,22,32)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_6_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)

		* Coef on event time 7 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(3,13,23,33)
		matselrc Vcov tempVcov, r(3,13,23,33) c(3,13,23,33) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2013+ $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_7_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 8 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(4,14,24,34)
		matselrc Vcov tempVcov, r(4,14,24,34) c(4,14,24,34)
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_8_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 9 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(5,15,25,35)
		matselrc Vcov tempVcov, r(5,15,25,35) c(5,15,25,35) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_9_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 10 for cohorts 2011, 2012, 2013, 2014
		* zero

		* Coef on event time 11 for cohorts 2011, 2012, 2013, 2014
		matselrc biw_long tempb, c(7,17,27,37)
		matselrc Vcov tempVcov, r(7,17,27,37) c(7,17,27,37) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_11_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 12 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,18,28,38)
		matselrc Vcov tempVcov, r(8,18,28,38) c(8,18,28,38) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2011 + $c_2012/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2012 + $c_2013/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2013 + $c_2014/($c_2011+$c_2012+$c_2013+$c_2014)*femx_12_2014
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 13 for cohorts 2011, 2012, 2013
		matselrc biw_long tempb, c(8,19,29)
		matselrc Vcov tempVcov, r(8,19,29) c(8,19,29) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012+$c_2013)*femx_13_2011 + $c_2012/($c_2011+$c_2012+$c_2013)*femx_13_2012 + $c_2013/($c_2011+$c_2012+$c_2013)*femx_13_2013
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 14 for cohorts 2011, 2012
		matselrc biw_long tempb, c(9,19)
		matselrc Vcov tempVcov, r(9,19) c(10,19) 
		matrix temp = tempb*tempVcov*tempb'	
		lincom $c_2011/($c_2011+$c_2012)*femx_14_2011 + $c_2012/($c_2011+$c_2012)*femx_14_2012 
		matrix b = b \ r(estimate)
		matrix V = V \ (r(se)^2 + temp)	

		* Coef on event time 15 for cohorts 2011
		lincom femx_15_2011
		matrix b = b \ r(estimate)
		matrix V = V \ r(se)^2		
		
		
		mat b = b[1..5,.] \ (0) \b[6..10,.]
		mat V = V[1..5,.] \ (0) \V[6..10,.]	
		mat b_2011=b_2011'
		mat b_2012=b_2012'
		mat b_2013=b_2013'
		mat b_2014=b_2014'	
		mat V_2011=V_2011'
		mat V_2012=V_2012'
		mat V_2013=V_2013'	
		mat V_2014=V_2014'		
		mat b_2011 = b_2011[1..5,.] \ (0) \b_2011[6..10,.]
		mat V_2011 = V_2011[1..5,.] \ (0) \V_2011[6..10,.]	
		mat b_2012 = b_2012[1..5,.] \ (0) \b_2012[6..10,.]
		mat V_2012 = V_2012[1..5,.] \ (0) \V_2012[6..10,.]	
		mat b_2013 = b_2013[1..5,.] \ (0) \b_2013[6..10,.]
		mat V_2013 = V_2013[1..5,.] \ (0) \V_2013[6..10,.]	
		mat b_2014 = b_2014[1..5,.] \ (0) \b_2014[6..10,.]
		mat V_2014 = V_2014[1..5,.] \ (0) \V_2014[6..10,.]	
		
		*mat biw_hrs = b,V,b_2011',V_2011',b_2012',V_2012',b_2013',V_2013'
		mat biw_hrs_m = b,V,b_2011,V_2011,b_2012,V_2012,b_2013,V_2013,b_2014,V_2014
		
		// output results to excel spreadsheet
		clear
		global n_obs = 11
		set obs $n_obs
		gen y = _n-5 // generate relative year

		* two-way FE estimates, young
		svmat b2way_hrs_m 

		gen ci1_u_m = sqrt(b2way_hrs_m2)*1.96+b2way_hrs_m1
		gen ci1_l_m = -sqrt(b2way_hrs_m2)*1.96+b2way_hrs_m1
		gen sd1_m = sqrt(b2way_hrs_m2)

		* IW estimates
		svmat biw_hrs_m

		gen ci2_u_m = sqrt(biw_hrs_m2)*1.96+biw_hrs_m1
		gen ci2_l_m = -sqrt(biw_hrs_m2)*1.96+biw_hrs_m1
		gen sd2_m = sqrt(biw_hrs_m2)
		* CATTs estimates are biw_hrs3 - biw_hrs7
		gen sd2011_m = sqrt(biw_hrs_m4)
		gen sd2012_m = sqrt(biw_hrs_m6)
		gen sd2013_m = sqrt(biw_hrs_m8)
		gen sd2014_m = sqrt(biw_hrs_m10)
		order y b2way_hrs_m1 biw_hrs_m1 biw_hrs_m3 biw_hrs_m5 biw_hrs_m7 biw_hrs_m9 sd1_m sd2_m sd2011_m sd2012_m sd2013_m sd2014_m
		
		export delimited using $data/weighted_mover3.csv, replace		
			
		appendfile $data/weighted_mover1.csv $data/weighted_mover2.csv 
		appendfile $data/weighted_mover2.csv $data/weighted_mover3.csv
		
	* Make graph
	clear all
	program varprocess
	{
		
		g l = _n-6

		keep if l<=5
		
		rename b2way_hrs1 FE
		rename biw_hrs1 IW
		rename biw_hrs3 CATT_1
		rename biw_hrs5 CATT_2
		rename biw_hrs7 CATT_3
		rename biw_hrs9 CATT_4	
		
		rename b2way_hrs_m1 FE_m
		rename biw_hrs_m1 IW_m
		rename biw_hrs_m3 CATT_1_m
		rename biw_hrs_m5 CATT_2_m
		rename biw_hrs_m7 CATT_3_m
		rename biw_hrs_m9 CATT_4_m
		
		rename sd1 FE_sd
		rename sd2 IW_sd
		rename sd2011 CATT_1_sd
		rename sd2012 CATT_2_sd
		rename sd2013 CATT_3_sd
		rename sd2014 CATT_4_sd
		
		rename sd1_m FE_m_sd
		rename sd2_m IW_m_sd
		rename sd2011_m CATT_1_m_sd
		rename sd2012_m CATT_2_m_sd
		rename sd2013_m CATT_3_m_sd
		rename sd2014_m CATT_4_m_sd
	
		rename b2way_hrs_old1 FE_old
		rename biw_hrs_old1 IW_old
		rename biw_hrs_old3 CATT_1_old
		rename biw_hrs_old5 CATT_2_old
		rename biw_hrs_old7 CATT_3_old
		rename biw_hrs_old9 CATT_4_old
	
		rename sd1_old FE_old_sd
		rename sd2_old IW_old_sd
		rename sd2011_old CATT_1_old_sd
		rename sd2012_old CATT_2_old_sd
		rename sd2013_old CATT_3_old_sd
		rename sd2014_old CATT_4_old_sd

		foreach var in FE IW FE_m IW_m CATT_1 CATT_2 CATT_3 CATT_4 CATT_1_m CATT_2_m CATT_3_m CATT_4_m FE_old IW_old CATT_1_old CATT_2_old CATT_3_old CATT_4_old {
			g upperCI_`var' = `var' + `var'_sd*1.96
			g lowerCI_`var' = `var' - `var'_sd*1.96
		}
	}
	end

	program scatterwCI 
		args var pat symb
	{
		
		global `var'_scatterwCI "rcap upperCI_`var' lowerCI_`var' l_`var', lw(medthick) lpattern(`pat') || scatter `var' l_`var' , c(l) msize(medlarge) m(`symb') "
	}
	end

	insheet using $data/weighted_mover3.csv, clear 

	varprocess
	gen l_IW = l - 0.05
	gen l_IW_old = l + 0.05
	gen l_IW_m	 = l
	scatterwCI IW dash T
	scatterwCI IW_old line O
	scatterwCI IW_m dash s

	twoway ${IW_scatterwCI} || ${IW_old_scatterwCI} || ${IW_m_scatterwCI} ///
		legend(order(2 "Never Move" 4 "Ever Moves" 6 "Moves Post")) xla(-5(1)5) ///
		xtitle("Years from CBA Extension") ytitle("Estimate: Log Salary")  ///
		graphregion(fcolor(white)) scheme(sj) 
	graph export "$out/AS_movers_weighted.png", replace


	
	
